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The success of the moving puncture method for the numerical simulation of black hole systems 
can be partially explained by the properties of stationary solutions of the 1+log coordinate condi- 
tion. We compute stationary 1+log slices of the Schwarzschild spacetime in isotropic coordinates 
in order to investigate the coordinate singularity that the numerical methods have to handle at the 
puncture. We present an alternative integration method to obtain isotropic coordinates that sim- 
plifies numerical integration and that gives direct access to a local expansion in the isotropic radius 
near the puncture. Numerical results have shown that certain quantities are well approximated by 
a function linear in the isotropic radius near the puncture, while here we show that in some cases 
the isotropic radius appears with an exponent that is close to but unequal to one. 
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I. INTRODUCTION 

The moving puncture method [U [2] is the basis of many 
of the recent successful simulations of black hole binaries 
in numerical relativity. For black hole puncture data, the 
term puncture refers to a single point on a hypersurface 
where the metric has a coordinate singularity that char- 
acterizes the presence of a black hole, while the physical 
singularity is not part of the hypersurface. 

The original puncture data ^31 HI [3 Ej uses the Brill- 
Lindquist "wormhole" topology [71|S], where the hyper- 
surface connects an outer with an inner asymptotically 
flat region, and the coordinate singularity occurs when 
compactifying the inner end by a singular conformal fac- 
tor in isotropic coordinates. For the Schwarzschild solu- 
tion in isotropic coordinates with radius r the conformal 
factor is i/) = 1 + 1^, where M is the mass, and the punc- 
ture singularity occurs for r ^ and Schwarzschild areal 
radius i? — > oo. The moving puncture method typically 
uses wormhole data as initial data, however the gauge 
conditions (1-l-log slicing and Gamma- freezing shift) lead 
after some gauge evolution to locally stationary solutions 
that are perhaps better described by the term "trum- 
pet" [HI Uni [m [H]- Trumpets for a stationary maxi- 
mal slice are discussed in [131 Ej • For Schwarzschild in 
isotropic coordinates, the stationary l-|-log slice is char- 
acterized by '0 ~ ^ for small r, and the trumpet ends 
at a finite Schwarschild radius, i? — > i?o ~ 1.312Af for 
r ^ 0. 

The analytic stationary l-l-log solution was first de- 
scribed for moving punctures in [3] (see [TS] for related 
results that however do not address the issue of black 
hole punctures). The analytic solution has been used to 
develop quite a complete picture of the geometry and 
global properties of the Schwarzschild solution for sta- 
tionary l-|-log slices [TT| [T2I [16], and detailed infor- 
mation about coordinate effects near the puncture has 
also been obtained through the analysis of ID numerical 
simulations in spherical symmetry |10L IllL 117] . 

In [9], we also discussed regularity at the puncture in 
terms of a Taylor expansion in the isotropic coordinate 



radius r around r = 0, similar in spirit to previous work 
on fixed punctures |181 [l9j . The motivation for such an 
investigation is to find out what sort of singularity a typi- 
cal 3D BSSN code in Cartesian coordinates has to handle 
at the puncture. On the one hand, the analytic results 
have shown that a solution for 1-1- log slicing exists, and it 
is an experimental fact that numerical codes are able to 
approximate this solution. On the other hand, it is not 
yet fully understood how the finite difference codes suc- 
ceed in obtaining accurate approximations to a black hole 
puncture that evidently is not regular at the puncture. 
For example, in the numerical coordinates the conformal 
factor and the lapse are not smooth at the puncture, and 
the BSSN extrinsic curvature displays a finite jump dis- 
continuity 'W . The analysis in [9J suggested that at least 
some of the leading order behavior near r = could be 
reliably obtained by a Taylor expansion since it approx- 
imates the numerical results rather well, although the 
expansion was at least partially an ad hoc ansatz. The 
existence of a solution with a particular small r behavior 
was assumed, and its consistency and some consequences 
were derived by inserting the ansatz into the full BSSN 
system and the gauge conditions. 

The purpose of the present paper is to derive the small 
r series in isotropic coordinates from the analytic solu- 
tion for stationary 1-1- log slices of the Schwarzschild solu- 
tion. The focus is on the calculations since as mentioned 
above the geometric picture has already been discussed 
elsewhere. We assume that isotropic coordinates approx- 
imate the 3D numerical simulations well near the punc- 
ture, although in fact the 3D numerical coordinates are 
not isotropic. We leave the comparison to actual sim- 
ulations, as well as the question how finite differencing 
works in this case to future work. The key question ad- 
dressed here is what the analytic results imply for the 
singularities at the puncture in isotropic coordinates. 

The main result is that the leading order behavior in 
some quantities is indeed described by terms linear in r 
as assumed in the previous ansatz. In particular, im- 
posing stationarity of the conformal factor in the given 
coordinates, dtip = 0, led to the conclusion that due to 
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the combined evolution and gauge equations ip ~ Rq jr 
for small r, where i?o is the Schwarzschild areal radius 
at the puncture point. Furthermore, the lapse vanishes 
at the puncture, a ~ 0, and the shift vector satisfies 
/3* ~ rn', where rt* = jr. These results also hold in 
our new analysis. However, the analysis of the present 
paper shows that non-leading order terms involve factors 
of r to some non-integer power. Making the point, the 
previous ansatz was 



(1) 



while as we show here the analytic stationary l-|-log slice 
implies for isotropic coordinates 
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As a second result, the methods we develop for the inte- 
gration of the isotropy condition are valid for all r (not 
just locally near r = 0), and they allow a somewhat 
simpler way to compute initial data for a Schwarzschild 
puncture in the moving puncture gauge than was avail- 
able previously |12j . 

In Sec. [n] we establish our notation for the 
Schwarzschild solution in 3+1 form. In Sec. |III[ we 
summarize results for stationary 1+log slicing of the 
Schwarzschild solution. Sec. IV describes a novel way to 



perform various integrations required to obtain isotropic 
coordinates. In Sec. |V|we comment on numerical meth- 
ods to compute the l-l-log data. Sec. VI concludes with 
a discussion. 

We use geometric units G = c = 1, and we set M = 1 
in calculations, although for clarity we give factors of M 
in some places. 



II. SCHWARZSCHILD SOLUTION IN 3+1 
FORM 

The basic equations for the Schwarzschild solution in a 
3+1 decomposition can be written in a number of ways. 
We choose to start with coordinates that encompass both 
Schwarzschild coordinates, which are convenient for solv- 
ing the Einstein equations, and spatially isotropic coordi- 
nates, which are typically used for puncture initial data. 



A. Time-independent, spherically symmetric 
metric 

Consider the metric in 3+1 form assuming spherical 
symmetry and time independence. We introduce coordi- 
nates i, r, 0, and for which the metric takes the form 



+ ^dr2 + i?2(d02_^sin2 



(3) 



where the coefficients only depend on r. The Arnowitt- 
Deser-Misner (ADM) variables are the 3-metric gij and 
the extrinsic curvature Kij, plus the lapse a and shift /9% 
e.g. [20 . For our choice of coordinates a = a(r), = 
/3(r)/(r), g„ = l//(r)2, g,, = R(r)\ K„ = k(r)lj(r)\ 
and Kq0 ~ l{r)R^{r). Furthermore, — gggs\T?9 and 
K4!4! — Kgg sin^ 9, and all other components vanish. The 
function /3 is the norm of P\ = P'' 1^'' / P = grrP^'P"' = 
/3i/3*. With the normal vector to constant r surfaces 
normalized such that giju^n^ = 1, the extrinsic curvature 
is expressed using two functions k(r) and l(r), Kij = 
kniUj + l{gij ~ n-irij), with trace 

K =k + 2l. (4) 

In Schwarzschild radial gauge, the areal radius equals 
the coordinate radius, 



i?(r) 



(5) 



which for vanishing shift, P — 0, leads to Schwarzschild 
coordinates. Isotropic coordinates are given by a confor- 
mal factor ip{r) such that 



f 



(6) 



where dr^ + r^dH.^ is the flat metric in spherical coordi- 
nates. Isotropic coordinates are therefore given by / and 
R that satisfy 



and the conformal factor is 



i?' 



r 



r 

R' 



(7) 



(8) 



B. Solution of the Einstein equations 

The Schwarzschild solution for our choice of coordi- 
nates can be written as 



fR', 



k 
I 



2M 

IT' 



i?'' 

P 

R' 



(9) 
(10) 

(11) 
(12) 



These relations can be obtained by coordinate transfor- 
mation from the standard form of the Schwarzschild so- 
lution, but it is also instructive to derive them by solving 
the ADM equations starting with (|3| . Given the two met- 
ric coefficients /(r) and R{r), the four remaining quan- 
tities for lapse, shift, and extrinsic curvature are deter- 
mined. In fact, only a depends directly on the choice 
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of /, and only a involves the coordinate r since the re- 
maining equations can be expressed in terms of i? with 
= df3{R)/dR. Note that 
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2M 
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which shows that we are using the Killing lapse and shift. 

As already mentioned, we obtain Schwarzschild coordi- 
nates for R{r) = r and /3 = 0, which implies = 1 — 2^, 
f = a, and fc = / = 0. Here we are looking for stationary 
solutions for 1-1- log slicing with a > 0, which is given 
in terms of one equation that fixes the lapse, and which 
leads to a shift and extrinsic curvature which are not 
identically zero. We fix the remaining freedom by requir- 
ing either the Schwarzschild radial gauge, i? = r, or the 
isotropic gauge, / = r/R. 



III. STATIONARY 1+LOG SLICING OF THE 
SCHWARZSCHILD SOLUTION 

We summarize several results on the integration of the 
stationary 1-l-log condition [9l [HJ [12j [TBI El] while adding 
some details useful for our purposes, for example a brief 
derivation of the integrated stationary l-|-log condition 
in the more general coordinates of Sec. Illj 



A. Integration of the stationary 1+log condition 



The 1-l-log slicing condition is 



(14) 



If the shift term is absent, which has also been used in 
some numerical simulations, then the stationary slice is 
a maximal slice, if = 0, which we will not discuss here 
(see [131 El)- Manifest stationarity dta = implies 



Cpa = 2aK, 
which in our coordinates becomes 

fPa' = 2{k + 2l)a. 



(15) 



(16) 



This equation looks rather harmless, but for an unfortu- 
nate choice of variables, say when written for the con- 
formal factor in isotropic coordinates, it becomes a sec- 
ond order ODE in the metric coefficients that cannot be 
solved explicitly. However, an easy way to proceed is to 
eliminate /, k, and I (but not (3) from (16) using the 
Schwarzschild solution ([9| - ( |T2| ). Assuming aj3/R' ^ 0, 
the l-|-log condition becomes 



R' 



(3 R 

This equation can be explicitly integrated. 



(17) 



(18) 



for some constant of integration C. With given in 
terms of a and i? by ( 13 ), the result is 



Ce" = - 1 



2M 



R\ 



(19) 



The above calculation generalizes immediately to various 
other lapse conditions, including maximal slicing, har- 
monic slicing, or the more general Bona-Masso slicings, 
although in the latter case the integration may not be 
possible explicitly. For maximal slicing, K = 0, and the 
term Ce" in ( 19 1 is replaced by C. For stationary har- 
monic slicing, £f3a = a^K, and the left-hand-side of ( 19 1 
is Ca^. 1-l-log slicing is special in that the term Ce" 
in ( 19 ) introduces a non-trivial dependence on a. In 
[T^ I16j. these different slicings are obtained as special 
cases of a more general slicing condition. 

The integrated 1+log equation ( [l9| does not exphcitly 
depend on r but only on a{r) and Ryr), i.e. it is the same 
equation whether we introduce the Schwarzschild radial 
gauge i? = r or not, and it is equally valid for isotropic 
gauge. Also, we did not have to specialize to i? = r to 
perform the integration. 



B. The critical point, determining C 

The constant of integration C is determined by requir- 
ing regularity of a' for < a < 1. Indeed, although 
integrating \\7\ when written in terms of (3 is straight- 
forward, we have to note a regularity issue for the right- 
hand-side. Since a'{R) — a'{r)/R'{r), substituting for (3 



in (17 1 gives (with M = 1) 

6 -I- 4i?(a2 - 1) 



a'{R) 



(2 + R{a^ ~2a-l))R 



(20) 



The stationary 1-l-log condition in this form is the start- 
ing point for its solution in [9) [21] . The right-hand-side 
of this equation can become singular if the denominator 
vanishes. We assume that R > and R'{r) > 0, which 
holds for Schwarzschild radial gauge (i? — r, R'(r) = 1) 
and which is natural in the context of isotropic coordi- 
nates. Regularity of a'(i?) requires that if the denomina- 
tor in (20 1 vanishes, then the numerator has to vanish, 
too. There are two solutions, one of which does fall into 
the interval < a < 1. The corresponding "critical" 
values a — ar and R—Rr are 



ac = \/l0-3 

Rr = 



0.162, 



1.541. 



(21) 
(22) 



The integrated H-log equation also has to hold at the 
critical point (ac,i?c), which implies 



C = -Rlt 



1 

128 
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(3 + V10)''e' 



1.554. (23) 



If we do not impose regularity at the critical point, then 
the integration constant of the 1-l-log equation can take 
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on different values, which leads to various other solutions 
(see the discussion in [TBI). 

We can evaluate a'{R) at the critical point by 
I'Hopital's rule. Taking derivatives of the numerator and 
denominator in ( 20 ) results in a quadratic equation for 



a' (Re) with the two solutions 



a'{Rc) 



-2± ^10 + 3^ 



0.607,-1.613. (24) 



Therefore, if we integrate (|20| as a first-order differential 
equation for a{R) starting at i? = Rc, there are two pos- 
sible initial tangents a'{Rc), which lead to two different 
solutions. The first solution with a'{Rc) > is the stan- 
dard stationary slice with a(R) monotonically increasing 
from to 1, while the second solution is decreasing with 
R. See Fig. [l] 



C. Properties of R{a) and a{R) 

We have seen that the stationary 1+log slicing con- 
dition can be integrated explicitly to obtain an implicit 
equation ( |19[ ) for a as a function of R, which however 
cannot be solved explicitly. On the other hand, the same 
equation can be read as a fourth-order polynomial equa- 
tion for R in terms of a. 



(a2-l)i?4 + 2i?3_Ce" = 0, 



(25) 



which can be solved explicitly, although there are four 
roots to consider. 

For C determined by regularity at the critical point, 
Mathematica finds four roots Rri, Rr2, Rri, and 
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Rri, two of which are real for < a < 1. The expressions 
are somewhat unwieldy so we do not show them here. By 
inspection, see Fig. [T] we can define 



R{a) — Rr4{a) for a < ac, 
R{a) — Rr2{oi) for a > ac- 



(26) 



The labeling of the roots and of the branches is arbitrary 
and only reflects Mathematica's choices in this regard. 
We pick the solution that corresponds to a ^ 1 for R 
oo, and that runs smoothly through the critical point. It 
is straightforward to see that R{a) is at a^- ^(ct) 
given in (26 1 is continuous, and so is its first derivative 



based on our discussion of a'(Rc) in the previous section. 
Taking further derivatives of a'{R) given in (20) should 
allow us to establish smoothness at Rc- 

Let us now turn to a(i?), which is implicitly defined 
by (251 or by i? = R{a). It is instructive to consider the 
two limiting cases a — > and i? — > oo. 

For a = 0, we obtain a finite value 

i?o = i?(0), a(i?o)=0. (27) 
Ro can be obtained directly from ( pS] ) by setting a — 0, 
Rti~2Rl + C^ 0, (28) 
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FIG. 1: The two real-valued roots of the integrated 1-l-log 
equation, with labeling introduced by Mathematica in one 
particular example (top panel). The standard stationary 
1+log slice corresponds to the curve that starts at a = 
and Ro ~ 1.312M, passes smoothly through « 0.162 and 
Rc ~ 1.541A/, and continues for a ^ 1 to ii ^ oo. The other 
smooth branch is ruled out by the boundary conditions. For 
the standard branch, the inverse radius, S — 1/R, shows a 
rather linear dependence on a, running from 5*0 = I/-R0 to 
(bottom panel). 



which has four solutions, still too unwieldy to show in 
closed form. For the case we consider, 



Ro = i?(0) w 1.312. 



(29) 



The other real branch gives Rq w 1.661, see also Fig. [T] 
The function a{R) can be Taylor-expanded around 

R = Rq, 



a(i?) = iR-Ro)a'iRQ) + -{R-Rofa"{Ro) + .. . , (30) 

where the leading constant vanishes by definition, 
a(i?o) = 0, the linear term is directly given by the sta- 
tionary l-|-log condition, (20 ), and further coefficients can 



be obtained by differentiating (20 1. The n-th derivative 



5 



of a{R) is obtained from (20 1 in terms R and lower or- A change of integration variables gives 



der derivatives of a(R), which can therefore be computed 
iteratively for n — 1, 2, .... The first two coefficients are 



a'(i?o) = 
a"{Ro) = - 



6-ARo 
(2 — Ro)Ro ' 
4 (-9RI + 29Rl 



27i?o + 6) 



(2 - i?o)3i?2 



(31) 



(32) 



Since R^ ^ and Rq ^ 2, all the derivatives are regular 
at Rq. Hence, a{R) possesses a regular Taylor series at 
-Ro- For example, 



a'{Ro) ~ 0.832. 



(33) 



Let us also consider the limit of large R. The Taylor 
series of a{R) in l/R begins with 



O 



1 

(34) 

See also Fig. fll Expansion ( 34 1 follows from ( 25 1 written 
as a = (1 - ^ + with e" 

argument assumes that a = l + 0(l/i?). 
assumption, Eqn. (|34l 



ferentiation of ( 25 1 



e for a — !■ 1. This 
Without this 
can be obtained by implicit dif- 
The expansion (34 1 shows that the 



constant of integration C of the l+log equation is not 
determined by the condition that a — 1 for i? — > cx3 (C 
only enters at fourth order) . This condition already holds 
due to asymptotic flatness of the Schwarzschild solution 
(§ - 



IV. ISOTROPIC COORDINATES 

Given R{a) and a{R) as the solution of the station- 
ary l-|-log condition, we turn to the isotropy condition, 
which relates the coordinate radius r to the areal radius 
R and the lapse a. In our coordinates, spatial isotropy is 
implied by / = r/R, For the stationary, spherically 
symmetric solution we have a — fR', ([9]), so the isotropy 
condition becomes 



a{r) 



rR'{r) 
R{r) ■ 



(35) 



This is one ODE involving two unknown functions a(r) 
and R{r), which however are related through a = a{R), 
or equivalently R — R{a) . 



A. Explicit integrals for r{R) and r(a) 
Written as 



dr _ dR 
~ ~ a(i?)i?' 



the isotropy condition leads to the integral 



(36) 



(37) 



r{a) = C2 exp 



1 dR{a) 
aR{a) da 



da. 



(38) 



The Ci are constants of integration. 



The task is to find practical methods to evaluate (37 1 
: (38). In contrast to maximal slicing [2], the inte- 



gration cannot be performed explicitly (i.e. Mathemat- 
ica does not know how, and the form of a{R) and R{a) 
makes the existence of an explicit solution unlikely) . 

is that the in- 



An immediate issue with (37 1 and (38 



tegral becomes divergent at its lower and upper bounds, 
a — Q and R — Rq, and a — 1 and R 00. In [l2l [T6]. 
we described a method that allows the accurate numer- 
ical evaluation of ( pS] ) in Mathematica. The integral is 
performed as which is split into two parts treating the 
cases a < 0.1 and a > 0.1 separately. For a > 0.1, the 
method is based on one partial integration t hat pulls out 
a factor of i?^/" on the right-hand-side of (38 1, so that 
the far limit r/R —f 1 for a — > 1 is directly implemented. 
The numerical integration relies on Mathematica's abil- 
ity to handle the singularity in the integrand for a — » 
automatically. 

Here we discuss a modified analytic formulation that 
regularizes the integral at both bounds by explicitly ex- 
tracting the problematic factors. This aids numerical 
evaluation because no special methods for divergent in- 
tegrals have to be employed, which we explore in Sec. \V\ 
Furthermore, the regularized expression allows an ana- 
lytic discussion of the limit a — > 0. Our previous method 
is not convenient for such an analysis because of the sin- 
gular limit and because the integration is split into two 
pieces. 



For large R, the integrand in (37 1 asymptotes to 1/i? 



1 



1 
R 



for R 



00, a 



(39) 



The integral J j^dR ~ InR is divergent, but in such a 
manner that exp f ^ dR — R and hence r ^ R. We 
therefore rewrite (36 1 as 



-dr — 



1 



obtaining for the integral 



i?exp 



1 -a 
aR 



dR, 



(40) 



(41) 



where we have introduced explicit integration limits, and 
we have fixed the constant of integration, Ci — 1. For 
large R we have the expansion (34 1 for a(i?), a ~ 1 — 



j^, and the integrand now has the asymptotic behavior 
~ . The integral is convergent for a fixed R > Rq, 

and 



R for i? ^ 00. 



(42) 
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For an alternative derivation of (41 1, note that a nat- 
ural variable to consider is r/R — f = ip~'^. Differ- 
entiating r/R with respect to R and using the isotropy 
condition, we get 



d{r/R) 
r/R 



aR 



-dR, 



(43) 



in analogy to ( 36 1 , which integrates directly to ( 41 1 . In 



other words, scaling r by i? avoids the detour of regu- 
larizing a singular integral whose exponential gives the 
required factor R for r R. 

For R approaching its lower limit, the integrand of the 
upper-limit regularized expression (41 1 and the original 
integrand both have the pole 



1 



aR 



1 



1 



ai{R — Rq)Rq 



for R Rq, a ^ 0, 



(44) 

where we have used the leading order linear term of the 
expansion of a at i?o, (30 ), with oi — a'{Ro) given in (31 ) 
and ( 33 ) . If we subtract this pole as we did with 1/i? for 
R oo, the lower limit becomes regular, but the upper 
limit picks up the singularity ln(i? — i?o)- Instead we 
write (with one factor Rq replaced by R) 



dr 

r 



1 



1 



1 



aR R ai{R-Ro)R 
dR dR 



dR 



R ai{R — Rf))R 
With/7,^di?=^ln(l-f ), 



(45) 



r(i?) = i? (1-1^) cxpI{R), 



/(i?) 



1 - a 



jRo \ dR 
R — Rq J R 



where 7 = l/(aii?o), or with (31) 



2-i?, 



7 = 







6-4i?o 



0.916. 



(46) 
(47) 

(48) 



Our final expression for r{R), 
m = (...). The goal 



( 46 1 , involves an integral 
was to remove all singu- 



lar terms from the integrand and the integral, and it is 
straightforward to see that I{R) ~ for i? ^ cx3 and that 
I{Ro) is finite, although I{Ro) is not explicitly available. 
The numerical result is that exp I{R) varies monotoni- 
cally between about 1.155 and 1, see Sec. |IVC| 

Given r{R), we can compute r(a) as r(i?(a)). Another 
possibility is to start with ( 38 1 and remove the singularity 
from the integrand as we did for (p7|), which results in 



r(a) = R{a)a'' exp J{a), 
R'{a) 



J{a) = 



(1-a) 



R{a) 



da 



(49) 
(50) 



Here R'{a) 

(poi 



l/a'{R), so that from the 1-1- log condition 



IT 
R 



(a) 



R{a^-2a~l) R' 



Q + AR{a^ - 1) 



R 



(0) 



i?0 



6-4i?, 







= 7- 
(51) 

Again, we find that the integrand of J{a) is regular, J{a) 
is finite, and exp J(q;) is of order 1. When comparing 
r{R), (|46|, and r(a), (|49]), for R Rq and a ^ 0, there 
is an additional factor of i? ~ i?o- 



B. ODE integration for a{r) 

The preceding section shows how to obtain explicit in- 
tegrals for r(a), although e.g. for initial data we want to 
compute the inverse relation a{r). Using the chain rule. 



da{r) _ da{R) dR{r) 
dr dR dr 



(52) 



The first factor on the right-hand-side is given by the 
(non-integrated) 1+log condition (20 1, the second factor 
by the isotropy condition ( 35 1 . This equation determines 
a{r) by an ODE, 



(53) 



da a R(a) 
dr r R'(a) ' 

where R' / R is given in terms oi R{a{r)) and a(r) in ( |51[ ). 
Previously, we integrated this equation as / dr/r obtain- 
ing r{a), cmp. (38). Incidentally, writing this as / da/a 
does not give us ayr) directly, since we do not have R{r) 
available and therefore cannot perform the integration 
directly. 

For R ^ Rq and a — s- 0, the isotropy condition in the 
form a = rR'{r)/R implies that r ^ or R' {r) 0. 
We consider r — > 0. In this limit, we can evaluate a'{r) 
by using I'Hopital's rule on the right-hand-side of ( [53| . 
Since - a'(0), we find that 



a'(0) = a'(0) 



6-4i?o 
2-i?n ' 



(54) 



compare 
a'(0) = 



(511 



A priori there are two possibilities. If 
then (54) is satisfied for all Rq ^ 2. If 
a'{Q) ^ 0, then we conclude that i?o = | ~ 1.333, which 
accidentally is rather close to Rq « 1.312 determined by 
the regularity condition at Re- 

Therefore, if we start integrating the combined 
isotropy and l-|-log condition (53 1 at r = 0, there is one 



distinguished case that is singled out by the assumption 
that a(0) = and a'(0) ^ 0, e.g. a(r) = air + a2r^ + ... 
with ai 7^ 0. This implies i?o = |- Only when inte- 
grating the ODE for r ^ 00 do we discover that these 
slices do not pass through {ac, Rc) and do not have the 
appropriate far limit. 

Eqn. (54) for a' (0) = is, however, consistent with the 
ansatz 



a{r) 



(55) 
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FIG. 2: The regularized integrand and its exponentiated 
integral exp/(S') ocurring in the integration formula for the 
isotropic radius r, plotted versus S — 1/R from 5* = to 
So = l/Ro ~ 0.762/M. The resulting effect in r{R) is a factor 
of order unity that varies monotonically from 1 at i? = cx3 to 
1.155 at i? = Ro. 



consistent with the calculation of r(a) leading to (49 1 



In other words, the leading order behavior of a{r) near 
r = can be directly obtained from the ODE ( 53 1 , if we 



use the information about global properties obtained for 
a{R), in particular that i?o 7^ I- 



The constant of proportionality in ( 55 1 is not deter- 



mined locally but requires knowledge of a global integral 
and depends on the critical constant C. Although we 



cannot invert (46 1 and (49 1 to obtain R{r) and a{r) ex- 
plicitly for all r, from (46 1 and (30) for r — + we obtain 



1 

— ( 

7 



a(r) ~ -e-^^^")/-^ 



r 

Ro 



1/7 



(56) 



C. Some quantitative results for isotropic 
coordinates 

Fig. [2] shows numerical results for the integral involved 
in the calculation of r{R) based on (46 1 and (67). The 
main feature is that exp/(i?) varies monotonically and 
almost linearly between 

exp/(i?o) « 1.155, exp/(oo) = l. (57) 

Considered for R ^ Ro, 

exp I{R) = exp /(i?o) [1 + (i? ~ Ro)I'{Rq) + ...], (58) 



1.0 



1.5 



2.0 2.5 
R/M 
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FIG. 3: The isotropic radius r{R) and its first derivative 
r'{R). As expected, r{R) is approximately linear for R vary- 
ing from 7?o to 00. Its first derivative has a weak pole singu- 
larity at Ro with exponent 7 — 1 —0.084. 



where /'(i?o) = -75^(1 + fi) with oi = a'(i?o) and 02 = 
a"{Ro). 

Since exp I{R) is of order unity, r(R) is approximated 

by 



i?(l 



Ro 
~R 



)T = R^-^{R-RQy 



(59) 



For comparison, the isotropic radius of wormhole, 
fixed-puncture data is given by i? = (1 -|- i^^r^, or 



= - [R-l±R^{R-2)^ 



(60) 



Apart from the symmetry corresponding to the worm- 
hole, there is some structural similarity to (59), i.e. one 
could define = \ and Row = 2. 

For R 00, we have r ~ R, while for R — s- Ro, 



iR) 



{R-Ror'''expI{Ro) 



1.181(i?-i?o 



,0.916 



(61) 



with r(i?o) — 0. Therefore, since 7 « 0.916 happens 
to be rather close to unity, and since exp/(i?o) is only 
slowly varying, we expect r{R) to be well approximated 
by the linear term {R — Ro)^ over the entire range of R. 

Fig. [3] shows the quantitative result for r{R), and also 
for r'{R). Since 7 < 1, the first derivative has a weak 
pole singularity at Ro, 



r'{R)^{R~Ro) 



7-1 



1 



R — i?( 



0.084 



(62) 



The same singularity occurs in r'(a), the first derivative 
of r{a) ^ a"', see (49). As a quantitative measure of the 
singularity, we plot (62) in Fig. |4j As another example, 
^/2,o.084 ^ inipiies a; « 1.1 X lO^^^^ 

A numerical code using isotropic coordinates typically 
does not compute r(a), but for example 



a{r) 



„1.091 



a'(r) 



„0.091 



a"{r) 



-0.909 



(63) 
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FIG. 4: Singular behavior of {R - i?o)^"^ ^ {R ~ i?o)"° °** 
for R in the range from Ro to Ro + W~'^M (top panel), and 



vanishing of r 



-1 + 1/7 



near r — (bottom panel). 



If a numerical code relied on a'(0) = at r = 0, then 
some inordinate amount of resolution would be required, 
see Fig. [4] However, typically even a"{r) is regularized 
properly in a standard BSSN moving puncture code due 



to the factor of -0 



multiplying the second deriva- 



tive of the lapse in the evolution equations for the extrin- 
sic curvature. The examples (63) show that the lapse is 
slightly "more regular" at r = than a{r) ^ r. 



V. NUMERICS 

For various applications of the stationary l-|-log slices 
there remains the task to find numerical solutions for 
R{a) and a{R), and to compute the integrals of the co- 
ordinate transformation to isotropic coordinates to ob- 
tain i?(r) and a{r). We comment on suitable numerical 
methods in this section. 



a simpler strategy is to work with a general root find- 
ing routine, for example the bracketed Newton method 
of [53]. This allows us to find both R{a) and a{R), as 
well as inverting other implicit relations, i.e. we typically 
need such an algorithm anyway. We save the work of im- 
plementing the explicit but lengthy formula produced by 
Mathematica. Newton's method typically requires only 3 
to 6 iterations to reach double precision, taking less time 
than evaluating Mathematica's formula (although there 
are more compact formulas available). 

To implement Newton's method, it is convenient to in- 
troduce S — 1/R, which transforms R G [Ro,oo) to the 
finite interval S G [5*0,0). The integrated H-log condi- 
tion gives 



F(5,a) 

dF 

~dS 
dF 

da 



= -1 + 2S -Ce^S^ (64) 
= 2-4Ce"S'3 (65) 

= 2a-Ce"S"^ (66) 



For example, to obtain S{a), fix a and iterate S by New- 
ton's method for F{S,a) = with derivative (65 1. The 



bracketed Newton method of [23] combines the standard 
Newton method with a bisection method that is used 
whenever the Newton step appears to fail, and which in 
particular avoids that the iteration leaves a given inter- 
val. To initialize the method, we have to localize the 
root approximately, see Fig. [5] Although S{a) and R{a) 
are monotonic, F{S,a) as a function of S is not. We 
can bracket the root by setting the starting interval to 
[0, S'c] if a € [ttc, 1] and to [S'cS'o] if a G [0,ac], where 
Sc — l/Rc- Mathematica's FindRoot routine applied to 
the problem requires similar care (i.e. it fails if the root 
is not bracketed correctly). We also found that Newton's 
method works rather well for R on the unbounded in- 
terval, for which we can obtain upper and lower bounds 
for the root from (34 1. Simple bisection does the job. 



too, giving on average the expected three digits for ten 
iterations. 

In practice, root finding works well and is applicable 
for both R{a) and a{R). It can be straightforwardly and 
efficiently implemented using a simple Newton method 
[53] (and also with Mathematica's FindRoot), the only 
issue being that we have to initialize differently for a < 
Ur and a > Ur. 



A. Numerical computation of R[ct) and oi{R) 



B. Numerical integration for isotropic coordinates 



An explicit form for R{a) is readily obtained based 
on the roots of a fourth order polynomial, see (26). Im- 



plemented in Mathematica, we may have to increase the 
default working precision for a close to 1, and we have 
to remove a small complex part of the result at round-off 
level. 

When working outside a package like Mathematica (for 
concreteness think of an implementation in C or C-l — h). 



Studying the near and far limit of the isotropy con- 
dition has left the integrals I{R) or J(oi) for numerical 
evaluation. Although the integrands are regular, at the 
boundaries they are defined by one-sided limits, which 
is sometimes called "regular improper" boundaries. In 
practice it is often sufficient to apply "open" integration 
formulas, for which the integrands are required close to 
but not at the boundaries. Accuracy can deteriorate close 
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FIG. 5: Solving the integrated stationary 1-fIog condition by 
root finding. Shown is F{S,a) as a function of 5 = 1/R for 
a = 0.0 to 1.0 in steps of 0.1 (top panel) and a zoom-in for 
a — ac ~ 0.10 to Qc -I- 0.10 in steps of 0.05 (bottom panel). 
The maximal value of 5* is So ~ 0.762, the critical value is 
Sc ^ 0.649. The solution S{a) is defined by F{S, a) = 0, and 
there may be one or two solutions for a given a. The root can 
be bracketed by considering [0, Sc] if a G [etc, 1] and [Sc, So] 
if a G [0,ac]. 



ator and denominator. For the purpose of integration, 
split the integral into two pieces as required to handle 
ac as an improper boundary. (If this is not done explic- 
itly, the Romberg integration routine may place integra- 
tion points randomly close to ac, leading to unexpectedly 
large errors there.) 

When computing initial data in isotropic coordinates, 
r{R) and r(a) actually still have to be inverted to obtain 
R{r) and a(r). This can be achieved by root finding 
similar to what was discussed in Sec. |V A| so we do not 
discuss this here in further detail. 

Integrating ( [53| as an ODE gives an alternative 
method to obtain a{r) directly. Examples for useful 
general purpose integrators are Mathematica's NDSolve 
or the Bulirsch-Stoer routine described in [53]. From 
the point of view of using a minimal number of numer- 
ical algorithms, the ODE integration can be traded for 
Romberg integration or whatever else is used to obtain 
J (a). Also, in this case we can avoid using a root finding 
routine if we implement R{a) by explicit root formulas. 
Integrate (53 1 for a{r) by some ODE algorithm using the 
explicit root for R(a) on the right-hand-side. Given a(r), 
compute R{r) = R{a{r)) using explicit roots. All other 
quantities are then directly available in terms of a(r) and 
i?(r). 

The issue with integrating an equation for a{r) is that 
we need a convenient starting point. Both r — and 
1/r = are improper but regular limits. A standard way 
to proceed is to factor the singularity and use Taylor 
expansions to integrate away from the boundary. For 
a{r) (as opposed to R{r)), we also have to handle the 
critical point at Vc- Note that Tc is not directly available, 
so we cannot start the integration there without, say, 
performing the explicit integrals discussed above. On 
the other hand, since accuracy near Tc is crucial for a 
continuous derivative of a(r), starting at Tc seems the 
most promising strategy. 

To end with a concrete suggestion, integrate 



to an improper boundary, and if this is an issue, Taylor 
expansions could be used. As already pointed out, if only 
a numerical answer is required, we can rely on numerical 
integration routines as available in Mathematica that ac- 
tually handle improper as well as some types of singular 
boundaries. In this case the transformation to regular 
boundaries is not always necessary, cmp. j12J. 

Such routines also handle integrating to oo, as required 
in I{R). As an alternative, we can change the integration 
variable to S ~ l/R, 



S'ir) = 



S{r)a{S{r)) 



(68) 



as an ODE for S{r), letting the automatic step size con- 
trol of the integrator handle the limits. Compute a{S) 
by Newton's method from (64 1. To avoid issues with 



starting at improper boundaries, start the integration at 



= 0.30345204271479997, S'(rc) = 4(V10 - 3), (69) 



/(i?) 



l/B. 



1 



75 



5*0 ~ S 



dS 



(67) 



This formula can be directly evaluated using e.g. 
Romberg integration 

The domain of integration in J (a) is finite, but in this 
case there is a potential issue at a^, where the regularity 
of R'{a) relies on the cancellation of roots in the numer- 



cmp. (22 1. That is, there is one magic number which we 



provide for the convenience of the reader based on the 
explicit integrals, = r{Rc) using Mathematica with 
20 digit accuracy. Given S{r), compute R{r) — l/S{r), 
a{r) = a{S{r)) and f{r) = = r/R{r). 

Independently of how the relation between r and R was 
obtained, all other quantities required to specify confor- 



mally flat initial data in ADM form follow from ( 10 H( 12 1 
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VI. DISCUSSION 

We have computed the stationary 1+log shccs of the 
Schwarzschild solution in isotropic coordinates. The 
computation goes beyond what was akeady known by 
providing alternative integration methods of the isotropy 
condition that simplify numerical integration, and by giv- 
ing direct access to local expansions in the isotropic ra- 
dius r near r = 0. 

Let us emphasize that the 3D numerical evolutions do 
not use isotropic coordinates. Specifically, even though 
"0 is chosen to obtain a metric with uniform determinant, 
the metric is not diagonal and the conformal metric com- 
ponents can easily deviate by 15% or 20% from 1. The 
reason is that although the initial data is conformally 
flat, there is some significant initial gauge evolution in 
which e.g. the shift evolves from identically zero to the 
stationary Gamma- freezing shift. During this time the 
coordinates are in motion, and the final transformation 
between the coordinate r and the areal R depends on de- 
tails of the shift condition and e.g. also the initial value of 
the lapse. For example, the size of the shift damping pa- 
rameter 77 directly influences the final gradients in R{r). 
See [ini m] for the rather involved procedure to trans- 
form the non-isotropic, numerical coordinates to some 
standard coordinates. 

For the analytic stationary l-l-log solution in isotropic 
coordinates, we showed that near the puncture factors of 
r as well as r^^'^ play a role. Concretely, i/)"^ ~ r and 
(3"^ ~ r, while a ~ r^ °^\ 

This was not detected in the numerical 3D data since 
we did not search for a small deviation from a 



we can integrate r'{R) towards larger R, but only spe- 
ciflc choices lead to the standard slice (compare [16]). 
For the standard slice, we determine Rq from (28 1 as a 



,1.000 



hke 



,1.091 



The main open question in this regard 



is whether and how much the numerical gauge (which 
as pointed out is not isotropic) affects these exponents. 
What are the deviations from isotropy, do they matter at 
the puncture? In the original ansatz, g^j = Sij + O(r^), 
implying that locally at the puncture the coordinates are 
isotropic. If this holds for arbitrary gauge to leading 
order, the analysis of isotropic coordinates leading to a ~ 
r^/''' applies also to the numerical simulations. 

Furthermore, if stationary 1-1- log slices in isotropic co- 
ordinates are chosen as initial data, then our analysis of 
the coordinate singularity at the puncture can answer the 
question what singularity the numerical evolution meth- 
ods have to handle, if the non-advected Gamma-freezing 
shift condition is used that preserves the isotropy. The 
advected Gamma-freezing shift condition leads to differ- 
ent radial coordinates, which in principle could also be 
analyzed along the lines discussed here for the isotropic 
case. (See [MJ [5S] for different types of advection in the 
shift condition using either dt or dt — P^di as time deriva- 
tive.) 

Based on our discus sion o n integrating the l-|-log equa- 
tion for a'{r) in Sec. IV B one point to make is that a 
local Taylor expansion at r = cannot take into account 
global boundary conditions or the regularity condition 
at ac- Given any starting value for Rq for r(i?o) = 0, 



function of C, which we set to its critical value in (23 1 
based on regularity at Re- If we choose a different Rq, 
we get a different slice. Expanding the BSSN equations 
and the gauge conditions near r = assuming a ~ r 
rather than a ^ r^l^ therefore may well lead to a consis- 
tent local solution, which however is incompatible with 
the global solution we are looking for. More generally, 
it would be relevant to determine which leading terms of 
the expansion are independent of global properties. For 
example, the conformal factor and the shift are linear in 
r by construction since — ]^ and /?'' without 
any approximation. The leading orders in (p| arise from 
R ~ Rq and ^{Rq) ~ — 1. Only the next to leading 

order terms involve non-linear r^/'*' terms. 

It is somewhat ironic that the Schwarzschild areal ra- 
dius (which is not used in 3D numerical simulations) re- 
sults in a regular Taylor expansion at R = Rq, while the 
transformation to isotropic coordinates relabels points 
such that instead of r the natural variable is r^^'^ . This 
happens because of the logarithmic nature of the com- 
bined l-|-log and isotropy conditions. The numerics does 
not even use isotropic coordinates. Are there globally 
regular coordinates that remain linear in radius near the 
puncture? 

Part of the motivation to analyze the local properties 
for r = in spherical symmetry is to prepare a similar 
study in axisymmetry. For axisymmetric black holes we 
do not have an analytic, stationary 1+log solution, but 
we can write down a Taylor expansion. Axisymmetric 
moving puncture data include the two important cases 
of a spinning black hole and of a black hole with linear 
momentum. Based on numerical experiments, the results 
for single non-moving punctures generalize to spinning 
and/or moving punctures, but we do not know yet how 
this is reflected in a stationary l-|-log slice discussion. 

The present paper can also be the starting point for 
the construction of initial data for multiple black holes 
in the moving puncture gauge since some simplifications 
and additional details about the analytic single black so- 
lution have been obtained. Standard methods superim- 
pose single black hole solutions to obtain an ansatz for 
the solution of the constraints. In [21 , the conformal 
thin-sandwich method is considered and the punctures 
are removed by excision at the apparent horizon. It is 
demonstrated that stationary H-log slices for binaries 
constructed with helical Killing vector boundary are not 
asymptotically flat (if the data are non-axisymmetric). 
However, it may still be possible to construct asymptot- 
ically flat initial data starting from single moving punc- 
tures using a different method. 
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